Summary:
The registration framework only supports images with sitkFloat32 and sitkFloat64 pixel types (use the SimpleITK Cast() function if your image's pixel type is something else).
Successful registration is highly dependent on initialization. In general you can:
There are many options for creating an instance of the registration framework, all of which are configured in SimpleITK via methods of the ImageRegistrationMethod class. This class encapsulates many of the components available in ITK for constructing a registration instance.
Currently, the available choices from the following groups of ITK components are:
The SimpleITK registration framework supports several optimizer types via the SetOptimizerAsX() methods, these include:
The SimpleITK registration framework supports several metric types via the SetMetricAsX() methods, these include:
The SimpleITK registration framework supports several interpolators via the SetInterpolator() method, which receives one of the following enumerations:
import SimpleITK as sitk
from downloaddata import fetch_data as fdata
%matplotlib notebook
import gui
import registration_gui as rgui
import matplotlib.pyplot as plt
import numpy as np
import os
OUTPUT_DIR = 'output'
import numpy as np
from ipywidgets import interact, fixed
We first read the images, specifying the pixel type that is required for registration (Float32 or Float64) and look at them. In this notebook we use a CT and MR image from the same patient. These are part of the training data from the Retrospective Image Registration Evaluation (RIRE) project.
#fetching COVID dicom data from directory
reader = sitk.ImageSeriesReader()
dicom_names = reader.GetGDCMSeriesFileNames('../../MIDRC-RICORD-1A/MIDRC-RICORD-1A-419639-000082/08-02-2002-CT CHEST WITHOUT CONTRAST-04614/601.000000-COR 3X3-86740/')
reader.SetFileNames(dicom_names)
fixed_image = reader.Execute()
dicom_names = reader.GetGDCMSeriesFileNames('../../MIDRC-RICORD-1A/MIDRC-RICORD-1A-419639-000082/08-02-2002-CT CHEST WITHOUT CONTRAST-04614/604.000000-COR 3X3-11320/')
reader.SetFileNames(dicom_names)
moving_image = reader.Execute()
#done
#Casting images from 32-bit integer to 32-bit float
fixed_image = sitk.Cast(fixed_image, sitk.sitkFloat32)
moving_image = sitk.Cast(moving_image, sitk.sitkFloat32)
#resampledImage = sitk.Resample(img2_original, img1)
gui.MultiImageDisplay(image_list = [fixed_image, moving_image],
title_list = ['Fixed Image', 'Moving Image'],
figure_size=(8,5));
#this image intensity is much clearer and effective for using the CenteredTransformInitializerFilter
intensity_profile_image = moving_image
fig = plt.figure()
plt.hist(sitk.GetArrayViewFromImage(intensity_profile_image).flatten(),bins=100);
intensity_profile_image = fixed_image
plt.hist(sitk.GetArrayViewFromImage(intensity_profile_image).flatten(),bins=100);
Estimate a 3D rigid transformation between images of different modalities.
We have made the following choices with respect to initialization and registration component settings:
We initialize registration by aligning the centers of the two volumes. To qualitatively evaluate the result we use a linked cursor approach, click on one image and the corresponding point is added to the other image.
#this allows us to map points to see
initial_transform = sitk.CenteredTransformInitializer(fixed_image,
moving_image,
sitk.Euler3DTransform(),
sitk.CenteredTransformInitializerFilter.GEOMETRY)
gui.RegistrationPointDataAquisition(fixed_image, moving_image, figure_size=(9,4), known_transformation=initial_transform);
Run the next cell three times:
registration_method = sitk.ImageRegistrationMethod()
# Similarity metric settings.
registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
registration_method.SetMetricSamplingPercentage(0.01)
#Interpolator
registration_method.SetInterpolator(sitk.sitkLinear)
# Optimizer settings.
registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100, convergenceMinimumValue=1e-6, convergenceWindowSize=10)
registration_method.SetOptimizerScalesFromPhysicalShift()
#uncomment line above, it makes all the differnce
# Don't optimize in-place, we would possibly like to run this cell multiple times.
registration_method.SetInitialTransform(initial_transform, inPlace=False)
# Connect all of the observers so that we can perform plotting during registration.
registration_method.AddCommand(sitk.sitkStartEvent, rgui.start_plot)
registration_method.AddCommand(sitk.sitkEndEvent, rgui.end_plot)
registration_method.AddCommand(sitk.sitkMultiResolutionIterationEvent, rgui.update_multires_iterations)
registration_method.AddCommand(sitk.sitkIterationEvent, lambda: rgui.plot_values(registration_method))
final_transform = registration_method.Execute(fixed_image, moving_image)
# Always check the reason optimization terminated.
print('Final metric value: {0}'.format(registration_method.GetMetricValue()))
print('Optimizer\'s stopping condition, {0}'.format(registration_method.GetOptimizerStopConditionDescription()))
Final metric value: -1.091190129503746 Optimizer's stopping condition, GradientDescentOptimizerv4Template: Convergence checker passed at iteration 54.
Qualitatively evaluate the result using a linked cursor approach (visual evaluation):
gui.RegistrationPointDataAquisition(fixed_image, moving_image, figure_size=(8,4), known_transformation=final_transform);
If we are satisfied with the results, save them to file.
moving_resampled = sitk.Resample(moving_image, fixed_image, final_transform, sitk.sitkLinear, 0.0, moving_image.GetPixelID())
sitk.WriteImage(moving_resampled, os.path.join(OUTPUT_DIR, 'RIRE_training_001_mr_T1_resampled.mha'))
sitk.WriteTransform(final_transform, os.path.join(OUTPUT_DIR, 'RIRE_training_001_CT_2_mr_T1.tfm'))
#fetching COVID dicom data from directory
reader = sitk.ImageSeriesReader()
dicom_names = reader.GetGDCMSeriesFileNames('../../MIDRC-RICORD-1A/MIDRC-RICORD-1A-419639-000082/08-02-2002-CT CHEST WITHOUT CONTRAST-04614/601.000000-COR 3X3-86740/')
reader.SetFileNames(dicom_names)
fixed_image = reader.Execute()
dicom_names = reader.GetGDCMSeriesFileNames('../../MIDRC-RICORD-1A/MIDRC-RICORD-1A-419639-000082/08-02-2002-CT CHEST WITHOUT CONTRAST-04614/604.000000-COR 3X3-11320/')
reader.SetFileNames(dicom_names)
moving_image = reader.Execute()
#cast imported images to floats
#Most often images have a high dynamic range. Thus, to write them to file in a format appropriate
#for use in a manuscript we will need to map the intensities to a low dynamic range (e.g. [0,255])
#In SimpleITK this is readily done with the IntensityWindowingImageFilter.
fixed_image = sitk.Cast(sitk.IntensityWindowing(fixed_image), sitk.sitkFloat32)
moving_image = sitk.Cast(sitk.IntensityWindowing(moving_image), sitk.sitkFloat32)
#changing the window level
ct_window_level = [835,230]
mr_window_level = [635,162]
gui.MultiImageDisplay(image_list = [fixed_image, moving_image],
title_list = ['fixed', 'moving'], figure_size=(8,4), window_level_list=[ct_window_level, mr_window_level]);